(N 



< 



(N 



c^ 



An efficient strategy to suppress epidemic explosion in 
heterogeneous metapopulation networks 

Chuansheng Shen^'^, Hanshuang Chen^, and Zhonghuai Houij 

^Hefei National Laboratory for Physical Sciences 

at Microscales & Department of Chemical Physics, 



O, 

CN , University of Science and Technology of China, Hefei, 230026, China 

^ ' '^Department of Physics, Anqing Teachers College, Anqing, 246011, China 



^School of Physics and Material Science, 
Anhui University, Hefei, 230039, China 



f. Abstract 

Q . We propose an efficient strategy to suppress epidemic explosion in heterogeneous metapopulation 
c/2 ■ 

f/3 ' networks, wherein each node represents a subpopulation with any number of individuals and is 

O ' 

*^ ' assigned a curing rate that is proportional to k" with k the node degree and a an adjustable 

>>: 

,S^ • parameter. We have performed stochastic simulations of the dynamical reaction-diffusion processes 
associated with the susceptible-infected-susceptible model in scale-free networks. We found that 

^ I the epidemic threshold reaches a maximum when the exponent a is tuned to be aopt — 1-3. This 

o\ ■ 

[^^ ■ nontrivial phenomenon is robust to the change of the network size and the average degree. In 
addition, we have carried out a mean field analysis to further validate our scheme, which also 

^—s I demonstrates that epidemic explosion follows different routes for a larger or less than Oopt- Our 



work suggests that in order to efficiently suppress epidemic spreading on heterogeneous complex 
networks, subpopulations with higher degrees should be allocated more resources than just being 



^ ', linearly dependent on the degree k. 
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I. INTRODUCTION 



In the last two decades, we have witnessed dramatic advances in complex networks re- 
search, which has been one of the most active topics in statistical physics and closely related 
disciplines [l|-l3|. The central issue in this field is to study how the network topology in- 



fiuences the dynamics 



Ml- 



As one of the typical dynamical processes built on complex 






networks, epidemic spreading has attracted more and more significant attention |7H24J|. De- 
spite much effort, many aspects of its role in the case of metapopulation models are still 
unclear and deserve more investigation. 

Very recently, metapopulation dynamics on heteregeneous networks, which incorporate 
mobility over the nodes, local interaction at the nodes, and a complex network structure, 
has gained great research attention 2J-|42|. In this context, reaction-diffusion (RD) pro- 
cesses have been widely used to model phenomena as diverse as epidemic and computer 



viruses spreading |24j-3C 



population evolution 



biological pattern formation 3lll32|], chemical reactions 33N35|. 
36| . and many other spatially distributed systems J37H40l| . In a series 
of important papers, Colizza et al. 25| provided an analysis of the basic RD process of 
the susceptible-infected-susceptible {SIS) model defined on heterogeneous metapopulation 
networks. Therein, each network node represents an urban area together with its population 
and edges represent air travel fiuxes along which individuals diffuse, coupling the epidemic 
spreading in different urban areas. They paid particular attention to the epidemic threshold 
Pc and found that pc is strongly affected by the topological fluctuations of the network for 
diffusing susceptible individuals. Later, Balcan and Vespignani J40[ extended such analy- 
sis to non-Markovian diffusive processes on complex networks, wherein individuals have a 
memory of their location of origin and displaced individuals return to their original subpop- 
ulation with a certain rate. Very recently, Vespignani \4:1\ reviewed and highlighted some of 
the recent progress in modelling dynamical processes that integrates the complex features 
and heterogeneities of real-world systems. Nevertheless, all the studies so far have treated 
the curing rate /x as a homogenous parameter, i.e., it is not dependent on the local property 
of the network node, such as the degree k. Note, however, in reality the curing rate of 
individuals should certainly be associated with the available medical resources in the local 
subpopulation, i.e., it is reasonable to assume that /i is a function of the degree k. It is 
therefore interesting to ask: how would the metapopulation dynamics of the SIS model, for 



instance, the epidemic threshold pc, depend on such a fc-dependent curing strategy? The 
answer to this question may provide useful instructions regarding the control of epidemic 
explosion in metapopulation networks. 

In the present paper, we have addressed such a question by considering a simple strategy, 
f^k ~ ^", where k denotes the node degree and a is an adjustable parameter. If a = 0, one 
recovers the usual cases studied in previous works. Herein, we mainly focus on the influence 
of varying a on pc- Interestingly, we found that pc bypasses a clear-cut maximum at a 
certain aopt, which corresponds to an optimal strategy to suppress epidemic explosion. This 
observation along with the value of aopt is robust to the change of the network size and the 
average network degree. To place the finding on a solid foundation, we have also performed a 
mean field (MF) analysis, wherein p^ is identified as the onset point where the global healthy 
state with no infected individuals loses stability. The MF equations successfully reproduce 
the Pc ~ « dependences, and also provide more insights regarding the routes to epidemic 
explosion for different values of a. 

II. MODEL DESCRIPTION 

We consider a system of A^ distinct subpopulations, each corresponding to a network 
node. Individuals inside each node run stochastically through the paradigmatic SIS model 



43l-l45| . Schematically, the stochastic infection dynamics is given by: 



S + 1A21, I^S (1) 

The first reaction reflects the fact that each susceptible individual becomes infected upon 
encountering one or more infected individuals at a probability rate /3. The second indicates 
that infected individuals are cured and become again susceptible at a fc-dependent rate 
Pk- Inside each network node, reaction processes take place under the assumption of a 
homogenous mixing and conserving the total number of individuals. After the reaction, 
individuals randomly diffuse along the edges departing from its local node. 

In this model, a significant and general result is that the system undergoes an absorbing- 
state phase transition with density p increasing, in analogy with critical phenomena |6|. 
Here p is defined as the total number of individuals divided by the number A^ of network 
nodes. The critical density pc indicates the epidemic threshold, what we are interested in. 



To begin, we perform our strategy on scale- free (SF) networks by using the Barabasi- 
Albert (BA) model [46| with power-law degree distribution p{k) ~ k~^. Scale-free networks 
are much more heterogeneous and serve as better candidates to test our strategy than other 
homogeneous networks, such as small-world or random networks. For a node i with degree 
ki, the curing rate is given by 

^'' ^ E~kf/N ^^^ 

Herein, /i^. is normalized such that the average curing rate remains constant: /i = -^ ^^ /ifc. = 
1. Note that in other related works about epidemic dynamics on networks, a fc-dependent 
strategy, but associated with the infection rate, had also been considered [1 71. l47l|. 



The system evolves in time according to the following rules J25| . The dynamics proceeds 
in parallel and considers a discrete time step representing the fixed time scale r of the 
process. The reaction and diffusion rates are therefore converted into probabilities. At 
each time step, the system is updated as follows. Inside each network node with degree 
k, each infected individual is cured and becomes an susceptible one with probability fik^. 
At the same time, each susceptible individual acquires infection from any infected one with 
probability 1 — (1 — /3t)"^, where nj is the total number of infected individuals in the 
node. After all nodes have been updated for the reactions, diffusion processes take place by 
allowing each individual to move into a randomly chosen neighboring node with probability 
DiT and Dst, for infected and susceptible individuals respectively, where Di{Ds) denotes 
the corresponding diffusion constant. In our simulation, the parameters are A^ = 1000, 
/3 = 0.5, Di = Ds = 1.0, T = 0.001 if not otherwise specified. Each plot is obtained via 
averaging over 20 independent simulation runs. 

III. SIMULATION RESULTS 

FigJTJ^a) shows how the proportion pi/p of infected individuals in the whole network 
increases with p, where pi = ^ rij/N denotes the density of infected individuals in the whole 
network, for four different values of a. Clearly, the system undergos a phase transition at 
a certain threshold density pc, above which pi/p monotonically increases from zero. For 
p < Pc the system stays in a 'healthy' state with pi = 0. Interestingly, pc reaches a largest 
value for a =1.3, compared to those for a =0, 1.0, 2.0. This is demonstrated more clearly 
in Fig{T] (b), where pc are plotted as a function of a for different network sizes N. The 



distinct peak locates at aopt — 1-3, which is rather robust to the change of network size N as 
shown in the inset. In addition, we have also investigated how this phenomenon depends on 
the average network degree {k). As shown in FigJT](c), the optimal value aopt also remains 
nearly constant with varying {k) from 4 to 14. 
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FIG. 1: (Color online) (a) The proportion pi/ p of infected individuals as a function of p for different 
a on 1000-node BA networks, (b) The epidemic threshold pc as a function of a for different network 
sizes N . The maximal threshold places Uopt — 1-3, which is indicated by vertical dotted line. The 
inset shows aopt as a function of N . All the networks have the fixed average network degree {k) = 6. 
(c)The epidemic threshold pc as a function of a for different {k). The inset shows aopt as a function 
of (/fc). iV = 4000. 



So far we have considered that all species diffuse with the same rate. In the following, we 
will take into account different diffusion rates for different species. For the sake of simplicity, 
we assume that infected individuals diffuse with a fixed rate Di = 1 and vary the diffusion 
rate of susceptible individuals Ds- The epidemic threshold pc as a function of a is plotted 
in Fig.2 for Ds = 0,0.005,0.05,0.5, and 1.0. Interestingly, the bell-shape dependence of Pc 
on a always exists for nonzero Ds, with the peak located at nearly the same optimal value 
ttopt — 1-3. The height of this peak decreases with Ds, and eventually pc is independent of 
a for Ds = 0. 

3.5 

3.0 

2.5 

Q." 2.0 
1.5 
1.0 
0.5 





1 






1 ' 1 




1 


□■ .. 


...□-.. 


■---Uy 


-a- 


'4^-->-^- 


- ■ o- ■ ■ 


■■■n 






9.' 


y 


/cOo ■ 


'"CX. 


■ 






/ (^ 




vcP^ '■ 


■«> ^ 












'O o ^' 


A 


"'X) 




.'O ^ 




1 


/ \ '^ 




■ 


0-' 


'.^' 




v/ 


\ ^' 




\ 


A'' 




t 


f 


■□■■■D3=0 V^ 
- o- D =0.005 \ 


\ 


A 


- 




yY 




s 


\^- 


- 




__ 


A 




-^- D =0.05 


/\ \ 






.V' 


/ 




s 




. 


^ 


' y 






-V- D =0.5 


N 


\ '•^-, 


" V" 


y> 






S 




^ " 


' <x^ 


^ 

1 






-O- 0^=1.0 

1 , 1 




1 



a 



FIG. 2: (Color online) The epidemic threshold pc as a function of a for different diffusion rates 
Ds- All the networks have the fixed (A;) = 6, A^ = 1000 and 7 = 3.0. 



Fig. 3(a) shows that p^ as a function of a for different infection rates /3. It can be found 
that the values of /3 do not influence the qualitative dependence of pc on a, i.e., a maximum 
Pc still shows up for the same optimal a. Nevertheless, the maximum pc corresponding to 
aopt do change with /3. In addition, we have also considered how the above findings depend 
on the network topology. To this end, we have performed simulations on SF networks 
with different exponents 7 and Erdos-Renyi (ER) random networks, pc as a function of a 
for different type of networks are shown in Fig.3(b). It is found that there still exists an 
optimal value of a, leading to the maximal threshold. For SF networks, the optimal value 
of a is always close to 1.3, while for ER networks, the pc ^ ex curve becomes not so sharp 
indicating that pc is not sensitive to the change of a. 
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FIG. 3: (Color online) The epidemic threshold pc as a function of a for different infection rates 
/3 (a) and for different network topologies (b). For (b), (3 = 0.5. All the networks have the fixed 
(A;) = 6 and Af = 1000. 

IV. MEAN FIELD ANALYSIS 

According to the stochastic simulation scheme, one may write down the following set of 
dynamic equations at a MP level, 



-Jr^ = Pi,k{f^Ps,k - Pk) + Di f k^p{k' \k)—pi^k' - Pi,k j 
—^ = PiAPk- /3ps,k) + Ds ( k^p{k'\k)—ps, 

\ k' 



k' — PS,k 



(3a) 
(3b) 



where pi^k and ps^k represent the average densities of infected and susceptible individuals, 
respectively, in the nodes with degree k. The first term in the right hand side of Eq.(3a) 



accounts for the change of infected individuals due to the reaction (infection and recovery) 
processes, and the second term accounts for the diffusion of infected individuals into and 
out of those nodes with degree k. Eq.(3b) can be interpreted in a similar manner. p{k' \k) 
represents the conditional probability that a node of degree k is connected to a node of 



degree k', which equals to k'p{k')/{k') 48|, |49| for BA networks. 



One notes that a thorough analysis of Eqs.(3) is not easy. For sake of simplicity, here 
we only consider the case Dj = Ds = 1. Then, substituting this into Eqs.(|3]) and using 
Pi = J2kPi^)Pi,k, one obtains 

-^ = Pi,k{l3ps,k - Pk) + j^pi - pi,k (4a) 

-^rr- = Pi,k{Pk - I^Ps.k) + -77TP5 - Ps.k (4b) 

For a = and thus /i^ = 1, it is already shown that pc = f Tpy [25]. But for o; 7^ 0, it 
is hard to get the explicit expression of pc from Eqs. (jl]) directly. Clearly, Eqs. (jl]) admit a 
steady state, which solves dpi^k/dt = dps,k/dt = 0, 



S,k 



which physically corresponds to the disease-free state. Intuitively, this healthy state will 
lose stability at the critical density pc, above which the steady state value of pi^x cannot be 
any more. Therefore, one can alternatively perform linear stability analysis of {p*.k,Ps.k) 
to get Pc- Following standard procedures, one can readily obtain the Jacobian matrix and 
calculate the eigenvalues {A}. The healthy state will lose stability when Xmax, the largest 
value of the real part of the eigenvalues, passes through zero from below. Note that explicit 
expression for Xmax is not available, but numerical calculation of it is easy. 

Fig. m (a) plots Xmax as a function of p for several values of a. The value of p where the 
Xmax = corresponds to pc- As expected, pc is the largest for a = 1.3 compared to those for 
other a. Fig. H] (b) presents pc as a function of a obtained from simulations (symbols) and 
MF analysis (solid line). Apparently, the MF results are in rather good agreements with the 
simulation ones in Fig.l. 

To get more insights into how the epidemic explosion takes place for different a, we turn 
to the eigenvector v = {{yi,k,^s,k)k=i,...} corresponding to Xmax at the onset of the phase 
transition, i.e., p = pc- The element v/ ^ of this vector measures the relative amplitude of 
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FIG. 4: (Color online) (a) and (b) correspond to the dependence of the largest eigenvalues \max 
on p for different a and the dependence of pc on a respectively, both on a synthesized 1000-node 
BA network and (A;) = 6. 

the fluctuation away from p} ^ = for nodes with given degree k. Therefore, the dependence 
of v/ fc on k quahtatively tells us how the epidemic explosion grows from the healthy state. 
In Fig. we depict the eigenvectors v/_fc as a function of k for different a. Interestingly, 
it can be observed that epidemic explosion starts from large-degree nodes for a less than 
ttopi, as shown by the "dotted", "dashed" and "dash dotted" lines in FigJSl while it is from 
small-degree nodes for a larger than acypt^ as shown by "solid" and "short dash dotted" lines. 
For a ~ c^opti "^i,k is not that sensitive on fc, indicating a relatively homogenous epidemic 
explosion. 

To reveal the underlying mechanism of the epidemic spreading for different a in more 
detail, we illustrate the time evolution of pi^k/ P-, the average density of infected individuals 




FIG. 5: (Color online) The eigenvectors v/ ^ corresponding to the dominant eigenvalue A^ 
as a function of k for different a. Other parameters are the same as in Fig. [H 



0, 



in the nodes with degree k in Fig. |6]for two particular values of a, one (a = 0.5) less than 
Uopt and the other (a = 2.5) larger than it. This can give us more detailed information 
about how the epidemic outbreak takes place on nodes with different degree k . We find 
that, for a = 0.5, the disease starts to spread from large-degree nodes, such as A;=88, 75 and 
60, as shown by the top three lines in Fig. 6(a); while for a = 2.5, the spreading starts from 
those nodes with relatively small degree, such as k=12, 17 and 8, as shown in Fig. 6(b). 
These phenomena indicate that there indeed exist two different epidemic explosion routes 
for a being less or larger than aopt, which are consistent with the analysis associated with 
the eigenvectors as shown in Fig. 5. 

The above different pathways regarding small or large a may be illustrated qualitatively 
in the following way. Consider the individuals in a given node are infected at the beginning. 
These patients will diffuse to neighboring nodes through the links. Certainly, nodes with 
larger degrees will have more chances to accept these patients. To efficiently suppress the 
epidemic explosion, the curing rates in such large-degree nodes should be relatively large to 
compensate these incoming patients via diffusion. Therefore, it is reasonable that fik should 
be an increasing function of k to maintain an effective epidemic control. Intuitively, one 
may imagine that the most efficient way is to keep linear dependence of fik on k, i.e., a = 1 
in our strategy, considering that every incoming patient via diffusion can be cured on time. 
However, this is not exactly the case because the reactions inside a node involve nonlinear 
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FIG. 6: (Color online) Time evolution pi^k/ P of the average density of infected individuals in the 
nodes with different degrees for BA model with N = 1000 and (A;) = 6. (a) a = 0.5, (b) a = 2.5. 

autocatalytic processes, which makes aopt larger than 1 (Unfortunately, why aopt is so robust 
to be about 1.3 is still open to us). If a is too large, which means that the medical resources 
are biased to large-degree nodes, the patients in small-degree nodes cannot be cured on time. 
In this case, disease will start to spread from those small-degree nodes. In the contrast case, 
the disease will start more abruptly from those large-degree nodes since the curing rates 
there are too smaller than required. These scenario are in agreement with the picture shown 
in Fig. |5] and [61 
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V. DISCUSSION AND CONCLUSIONS 

One should note that the a value can not be arbitrary for the real world, if we accept the 
power-law dependence. Following the recipe of Eq.2, for a scale free network with minimum, 
mean and maximum degree respectively of 2, 5 and 100, the recovery rate will range from 
0.4 to 20 in the simplest case of linear dependency [a = 1). This large difference is to some 
extent not reasonable, which implies that the optimal control with a = 1.3 is hard to be 
realized practically. Nevertheless, as a model study, we can just change a as we want to see 
what we can find. If, for instance, we tune a to a reasonable non-zero value, say a = 0.5, the 
ratio of the maximal and minimal /i would be about 10 for a network with k ranging from 
1 to 100, which can also lead to a much better epidemic control (pc = 1-37) than previous 
case of q; = {pc = 0.95). Therefore, our work has indeed provided an efficient strategy to 
suppress the epidemic explosion. 

In summary, we have studied a variant of SIS model defined on scale-free metapopulation 
networks, wherein the curing rate in a node with degree k is proportional to A;". By detailed 
numerical simulations, we show that the epidemic threshold reaches a maximum value when 
a is tuned to be aopt — 1-3, which corresponds to an optimal control strategy to suppress 
epidemic explosion and is robust to the change of network size or average degree. We have 
also performed a mean field analysis to further elucidate this strategy and unravel the distinct 
pathways to epidemic spreading for a larger or less than aopt- Out findings suggest that a 
proper allocation of medical resources can best suppress the epidemic explosion, which could 
be of great importance in practical epidemic control. 
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